Hydrological, Geomorphological, and Biological Applications

The topography of a catchment has a major impact on the hydrological, geomorphological. and biological processes active in the landscape. The spatial distribution of topographic attributes can often be used as an indirect measure of the spatial variability of these processes and allows them to be mapped using relatively simple techniques. Many geographic information systems are being developed that store topographic information as the primary data for analysing water resource and biological problems. Furthermore, topography can be used to develop more physically realistic structures for hydrologic and water quality models that directly account for the impact of topography on the hydrology. Digital elevation models are the primary data used in the analysis of catchment topography. We describe elevation data sources, digital elevation model structures, and the analysis of digital elevation data for hydrological, geomorphological, and biological applications. Some hydrologic models that make use of digital representations of topography are also considered. K E Y WORDS Basin topography Digital elevation models Terrain analysis Hydrologic models INTRODUCTION The demands placed on hydrologic models have increased considerably in recent times. Although it was once sufficient to model catchment outflow, it is now necessary to estimate distributed surface and subsurface flow characteristics, such as flow depth and flow velocity. These flow characteristics are the driving mechanisms for sediment and nutrient transport in landscapes and unless they can be predicted reasonably well, water quality models cannot be expected to adequately simulate sediment and nutrient transport. The control exerted by topography on the movement of water within the landscape is fundamental t o the prediction of these flow characteristics. A great deficiency of many hydrologic and water quality models currently in use is their inability to represent the effects of three-dimensional terrain on flow processes and the spatial variability of hydrologic processes without gross, and often unrealistic, simplifications. In response to this challenge, methods are being developed to digitally represent terrain in hydrologic models. In addition to these complex models, used largely for research into catchment processes, there is a demand for simpler techniques to assist with day-to-day land management. Action agencies responsible for land and water management in many parts of the world are being required to identify those areas of land susceptible to various types of environmental hazard and degradation such a s erosion, sedimentation, salinization, nonpoint source pollution, and water logging (Herndon, 1987) and to assess and manage biological productivity 0885 -6087/9 1/01 0003 -28$14.OO ( 1991 by John Wiley & Sons, Ltd. Received and accepted 6 September 1990 4 1. D. MOORE. R. B. GRAYSON AND A. R . LADSON and diversity within landscapes. A number of hydrologically-based, topographically-derived indices appear to be particularly powerful and useful for determining this susceptibility to hazard (Moore and Nieber, 1989). Many geographic information systems and resource inventory systems are being developed storing topographic information as primary data for use in analysing water resource and biological problems. This paper reviews the availability of digital elevation data and its accuracy, digital representation of topography, analysis of the digital data for hydrological, geomorphological, and biological applications, and describes some models that make use of digital representations of topography. We use examples from Australia and the United States of America in our discussions because of our familiarity with hydrologic research and practice in these two countries. DIGITAL ELEVATION MODELS A Digital Elevation Model (DEM) is an ordered array of numbers that represents the spatial distribution of elevations above some arbitrary datum in a landscape. I t may consist of elevations sampled at discrete points or the average elevation over a specified segment of the landscape, although in most cases i t is the former. DEMs are a subset of Digital Terrain Models (DTMs) which can be defined as ordered arrays of numbers that represent the spatial distribution of terrain attributes. The acquisition. storage, and presentation of topographic information is an area of active research. Generally, raw elevation data in the form of stereo photographs or field surveys, and the equipment to process these data, are not readily available to potential end users of a DEM. Therefore, most users must rely on published topographic maps or DEMs produced by government agencies such as the United States Geological Survey (U.S.G.S.) or the Australian Surveying and Land Information Group (AUSLIG-formerly the Division of National Mapping). This section describes the three major methods of structuring a DEM, the sources of elevation data in Australia and the United States, how DEMs are commonly produced and gives some indication of their quality. The section concludes with a discussion of some simple techniques for analysing elevation data for the estimation of topographic attributes such as slope and aspect and gives examples of methods that can be applied to each of the three DEM structures. Dutu networks ,for Digital Elevation Models When discussing the use of DEMs it is important to consider the way in which the surface representation is to be used. The ideal structure for a DEM may be different if it is used as a structure for a dynamic. hydrologic model than if it is used to determine the topographic attributes of the landscape. There are three principal ways of structuring a network of elevation data for its acquisition and analysis, as illustrated in Figure 1 . Triangulated Irregular Networks (TINS) usually sample surface specific-points, such Figure I . Methods of structuring an elevation data network: (a) square-grid network showing a moving 3 x 3 submatrix ccntrcd on node 5; (b) triangular irregular network-TIN; and (c) contour-bascd network DIGITAL TERRAIN MODELLING I : REVIEW OF APPLICATIONS 5 as peaks, ridges, and breaks in slope, and form an irregular network of points stored as a set of x, y, and z coordinates together with pointers to their neighbours in the net (Peucker et al., 1978; Mark, 1975). The elemental area is the plane joining three adjacent points in the network and is known as a facet. Grid-based methods may use a regularly-spaced triangular, square, or rectangular mesh or a regular angular grid, such as the 3 arc-second spacing used by the US. Defence Mapping Agency. The choice of grid-based method is related primarily to the scale of the area to be examined. The data can be stored in a variety of ways, but the most efficient is as z coordinates corresponding to sequential points along a profile with the starting point and grid spacing also specified. The elemental area is the cell bounded by three or four adjacent grid-points for regular triangular and rectangular grid-networks, respectively. Contour-based methods consist of digitized contour lines and are stored as Digital Line Graphs (DLGs) in the form of x, y coordinate pairs along each contour line of specified elevation. These can be used to subdivide an area into irregular polygons bounded by adjacent contour lines and adjacent streamlines (Moore, 1988; Moore and Grayson, 1989, 1990) and are based on the stream path analogy first proposed by Onstad and Brakensiek ( 1 968). The most widely used data structures consist of square-grid networks because of their ease of computer implementation and computational efficiency (Collins and Moon, 1981). However, they do have several disadvantages, including: (1) they cannot easily handle abrupt changes in elevation, although Franke and Nielson (1983) discuss possible ways of modelling these discontinuities; (2) the size of grid mesh affects the results obtained and the computational efficiency (Panuska et al., 1990); (3) the computed upslope flow paths used in hydrologic analyses tend to zig-zag and therefore are somewhat unrealistic; and (4) precision is lacking in the definition of specific catchment areas. Since regular grids must be adjusted to the roughest terrain, redundancy can be significant in sections with smooth terrain (Peucker et al., 1978), whereas triangulated irregular networks are more efficient and flexible in such circumstances. Olender (1980) compares triangulated irregular networks with regular grid structures. Topographic attributes such as slope, specific catchment area, aspect, plan, and profile curvature can be derived from all three types of DEMs. Methods of doing this and potential uses of these and other compound attributes are discussed in detail later. However, the most efficient DEM structure for the estimation of these attributes is generally the grid-based method. Contour-based methods require an order of magnitude more data storage and do not provide any computational advantages. With TIN structures, there can be difficulties in determining the upslope connection of a facet, although these can be overcome by visual inspection and manual manipulation of the TIN network (Palacios and Cuevas, 1989). The irregularity of the TIN makes the computation of attributes more difficult than for the grid-based methods. For dynamic hydrologic modelling, there are quite different considerations. Mark (1978) noted that grid structures for spatially partitioning topographic data are not appropriate for many geomorphological and hydrological applications. He stated that 'the chief source of this structure should be the phenomena in question, and not problems, data, or machine considerations, as is often the case'. Hydrologic models simulate the flow of water across a surface so the elemental areas of the DEMs should reflect this requirement. Contour-based methods have important advantages in this regard (Moore, 1988; Moore and Grayson, 1989, 1990) because the structure of their elemental areas is based on the way in which water flows on the land surface. Orthogonals to the contours are streamlines so the equations describing the flow of water can be reduced to a series of coupled one-dimensional equations. Digital elevation data sources Digital elevation data are available for selected areas in the United States from the National Cartographic Information Center (NCIC) in several forms (Dept. Interior-U.S.G.S., 1987). DEMs are available on a 30 m square grid for 7-5 minute quadrangle coverage, which is equivalent to the 1:24000-scale map series quadrangle. These data are produced by the U.S.G.S. for the NCIC from: ( I ) Gestalt Photo Mapper I1 (Kelly et al., 1977); (2) manual profiling from photogrammetric stereomodels; (3) stereomodel digitizing of contours; and/or (4) digital line graph data. DEMs based on a 3 arc-second spacing have been developed by the Defense Mapping Agency (DMA) for 1 degree coverage, which is equivalent to the 1 :250000-scale map series quadrangle. Line map data in digital form, known as Digital Line Graph (DLG) data, are available for 7.5 minute and 15 minute topographic quadrangles, the 1: 100000-scale quadrangle series, and 1 :20000006 1. D. MOORE, R. B. GRAYSON A N D A. R. LADSON scale maps. A 30 arc-second DEM is now available for the conterminous United States from the National Geophysical Data Center (NGDC) of the National Oceanic and Atmospheric Administration (NOAA). Moore and Simpson (1982) developed a coarse resolution DEM of Australia with a 6 minute (10 km) grid spacing and the Australian Survey and Land Information Group (AUSLIG) is currently developing a 1 : 1000000-scale grid DEM using a 18 arc-second grid-spacing (i.e. approximately 500 m) (Trezise and Hutchinson, 1986). Recently, Hutchinson and Dowling (1991) developed a 1.5 minute grid DEM of Australia. Both the 1.5 minute and 18 arc-second DEMs were developed using the procedures described by Hutchinson (1989a) which use both spot height data and stream line data. At the global scale, a 5 minute global DEM is available from NGDC-NOAA. As noted above, direct photogrammetric techniques are available to produce DEMs from stereophoto pairs. With the Gestalt Photomapping system, photo coordinates of control points are measured and their locations and elevations input (Kelly et al., 1977). A least squares solution is used to produce a stereo model that can be processed off-line to generate a DEM. Stereo images available from the French earth observation satellite, SPOT (Le Systeme Pour I’Observation de la Terre), can be used to produce orthophotos and DEMs in much the same way as conventional air photography. Satellite data has the advantage that it can be purchased in digital form and directly accessed by computers. I f digital elevation data for a particular study area are not available they can be derived by digitizing existing topographic maps. Contour lines can be digitized automatically using the processes of raster scanning and vectorization (Leberl and Olson, 1982) or by manually using a flat-bed digitizer and software packages. A variety of low-cost flat-bed digitizers and software packages are now available for most personal computers. Individual contours can be digitized and retained in contour form as DLGs or interpolated onto a regular grid or TIN. Similarly, spot heights can be digitized and analysed as an irregular network or interpolated onto a regular grid. Hutchinson (1988, 1989a) describes a new efficient finite difference method of interpolating grid DEMs from contour line data or scattered surface-specific point elevation data. The method is innovative in that it has a drainage enforcement algorithm that automatically removes spurious sinks or pits and calculates stream lines and ridge lines from points of locally maximum curvature on contour lines. The most common method of performing a contour to TIN transformation is via Delaunay triangulation (McLain, 1976), but this method can produce poorly configured triangular facets where a facet edge crosses a contour segment or all three vertices are on the same contour. Christensen (1987) developed a parallel pairs procedure as part of a medial axis transformation method that overcomes these problems. Digital elevution model quality Many published DEMs are derived from topographic maps so their accuracy can never be greater than the original source of the data. For example, the most accurate DEMs produced by the U.S.G.S. are generated by linear interpolation of digitized contour maps and have a maximum root mean square error (RMSE) of onehalf contour interval and an absolute error no greater than two contour intervals in magnitude (Dept. Interior U.S.G.S., 1987). Care must be exercised when fitting mathematical surfaces to DEMs as the resulting models give the appearance of ‘generating’ data, but the accuracy of the data is unknown. The U.S.G.S. DEMs are referenced to ‘true’ elevations from published maps that include points on contour lines. bench marks, or spot elevations (Dept. Interior-U.S.G.S., 1987). However, these ‘true’ elevations also contain errors. I t is therefore apparent that all DEMs have inherent inaccuracies not only in their ability to represent a surface but also in their constituent data and it ‘behoves users to become aware of the nature and the types of errors’ (Carter, 1989). The quality of the U.S.G.S. DEM data are classified as Level I , 2, or 3 (Dept. Interior-U.S.G.S., 1987). Level 1 is the standard format and has a maximum absolute vertical error of 50 m and a maximum relative error of 21 m. Most 7.5 minute quadrangle coverage is classified as Level 1. Level 2 data have been smoothed and edited for errors. DEMs derived by contour digitizing are classified as Level 2. These da!a have a maximum error of two contour intervals and a maximum root mean square error (RMSE) of one half of a contour interval. Level 3 data have a maximum error of one contour interval and a maximum RMSE of one third of a contour interval-not to exceed 7 m. The U.S.G.S. does not currently produce Level 3 elevation DIGITAL TERRAIN MODELLING I: REVIEW OF APPLICATIONS 7 data. The vertical accuracy of the Gestalt Photomapping system is 0.022 per cent to 0-03 per cent of the aircraft flying height and is therefore capable of producing a DEM of level 3 accuracy (Kelly et al., 1977). Konecny et a/. (1987) showed that SPOT stereoscopic data may be used for topographic mapping up to a scale of I :25000 and is capable of producing a DEM with a grid spacing of 500 m and a root mean square error of 7.5 m. The errors contained in the U.S.G.S. grid-based DEMs were classified by Carter (1989) as either global or relative. Global errors are systematic errors in the DEM that can usually be corrected by applying transformations including linear or nonlinear translation, rotation, and scaling to the whole DEM. Like Carter, the authors have found mismatching elevations along the boundaries of adjacent 7.5 minute DEMs to be the most pervasive global errors. Relative errors occur when a few elevations are in ‘obvious error relative to the neighbouring elevations’ (Carter, 1989). The errors in DLG data include artificial peaks and blocks of terrain extending above the surrounding terrain, particularly along ridge lines. These errors may be due to deficiencies in the interpolation algorithms that produced the DLG data or to an excessive spacing between sampled coordinate pairs along the line data. Thapa (1990) describes a method of optimal compression of DLG data that may overcome some of these problems. Analysis of elevation data The topographic attributes of a landscape can be calculated directly from a DEM using only the point values without the assistance of surface fitting, smoothing operations or the assumptions of continuity (Collins, 1975). This approach has limited usefulness, is restricted to grid-based DEMs, and does not produce physically realistic results particularly in the calculation of flow directions in flat areas (Douglas, 1986). The most common method of estimating topographic attributes involves fitting a surface to the point elevation data using either linear or nonlinear interpolation. The sequential steepest slope algorithm developed by Leberl and Olson (1982) is an example of a linear interpolation technique for generating grid elevation data. The interpolation is carried out firstly along ridge and drainage lines and then elevations are interpolated for all other points in the grid using ridge, drainage, and contour lines as lines of known elevation. Barnhill (1983) classifies nonlinear surface fitting methods in t w o ways: as local or global and patch or point schemes. Global methods utilize all or most of the elevation data to characterize the surface at a point and have the advantage of preserving continuity. Their chief disadvantage is their high computational cost which is proportional to n3, where n is the number of data points (Hutchinson, 1989a). For local methods, the fitted surface (often a simple polynomial) at a point depends only on nearby data and the computational costs are proportional to n (Hutchinson, 1989a). Patch surfaces consist of small curved patches that are joined together smoothly, whereas point methods construct the surface using only information given at discrete points. A wide variety of methods are available for fitting elevation surfaces to point data and some are described by Hutchinson (1984). They include: kriging (Gandin, 1963; Matheron, 1973; Delfiner and Delhomme, 1975; Jupp and Adomeit, 1981; and Dubrule, 1984), for which numerous commercial packages are available; local interpolation methods (surveyed by Barnhill and Boehm, 1983); moving average methods; and spline interpolation (Dubrule, 1984; Hutchinson, 1984). More recently, Hutchinson (1988,1989a) has developed an iterative finite difference interpolation method for use with irregularly distributed data that has the efficiency of a local method without sacrificing the advantages of the global methods. This method uses a nested grid strategy that calculates grid DEMs at successively finer resolutions. Most digital terrain analysis methods are based on grided data structures and for these local interpolation methods are the simplest and easiest to implement. This approach has been used by Evans (1980), Zaslavsky and Sinai ( 1 98 1 ), Mark ( I 983), Jenson (1985,1987), Zevenbergen and Thorne ( 1987), Jenson and Domingue ( 1 988), and Moore and Nieber (1989). One simple grid-based local interpolation method that appears very promising is Snyder et al.’s (1984) sliding polynomials. I t has the advantage over methods like Zevenbergen and Thorne’s (1987) in that the continuity of the fitted function and its first and second derivatives are preserved between adjacent patches. One problem with the analysis of digital elevation data for hydrologic applicatiow is the definition of drainage paths when the DEM contains depressions or flat areas. Some depressions are data errors while 8 1. D. MOORE, R . B. GRAYSON AND A. R . LADSON others are natural features or excavations (Jenson and Domingue, 1988; Hutchinson, 1989a). The hydrologic significance of depressions depends on the type of landscape represented by the DEM. In some areas, such as the prairie pot-hole region in the upper Midwest of the United States, surface depressions dominate the hydrologic response of the landscape. In areas with coordinated drainage, depressions are an artifact of the sampling and generation schemes used to produce the DEM. In landscapes with natural depressions, the numerical filling of depressions in the DEM is used as a method of determining storage volumes and to assign flow directions that approximate those occurring in the natural landscape once the depressions are filled by rainfall and runoff (e.g. Moore and Larson, 1979). Mark (1983) discusses the smoothing of data sets to remove depressions, but this has the disadvantage of affecting all the elevations in the data set and will not remove large depressions. O'Callaghan and Mark (1984) and Jenson (1987) proposed algorithms to produce depressionless DEMs from regularly spaced grid elevation data. If the depressions are hydrologically significant then their volume can be calculated. Jenson and Domingue (1988) used the depressionless DEM as a first step in assigning flow directions. Their procedure is based on the hydrologically realistic algorithm discussed by Mark ( 1 983) and O'Callaghan and Mark (1984), but is capable of determining flow paths iteratively where there is more than one possible receiving cell and where flow must be routed across flat areas. Hutchinson's (1988, 1989a) method, which was briefly described earlier, includes an automatic drainage enforcement algorithm that removes spurious sinks or pits. This can greatly simplify the task of obtaining a depressionless DEM. If a surface defined by the function F(x, y , z) is fitted to the DEM, then a number of hydrologically important topographic attributes can be derived from this function at the point (.yo, yo, zo). As an example. we will now demonstrate how two of these attributes, aspect and the maximum slope can be determined from the function F(x. y, z ) using elementary geometry. We will then briefly describe three elementary local-point methods of representing surfaces using triangulated irregular networks, square-grid networks, and contourbased networks, respectively. Thefitted s u t j i x ~ ~ The equation of the tangent plane to the point (xo. yo, zo) is